Esse notebook busca introduzir o algoritmo de Árvore de Decisão, para isso será usado o dataset Iris.
# estabelecendo o ambiente
suppressMessages({
library(tidyverse)
library(magrittr)
library(knitr)
library(GGally)
library(rpart)
library(rattle)
library(rpart.utils)
})
setwd("~/Dropbox/kaggle/iris-species/")
opts_chunk$set(cache=TRUE)
data(iris)
iris %<>% as_tibble()
O dataset contém medidas de comprimento e largura da pétala e da sépala de três espécies da flor Íris:
VersicolorSetosa
Virgínica
Entendendo a flor:
Uma amostra do dataset, ao todo ele contém 150 amostras, 50 de cada espécie.
iris %>% head() %>% print()
## # A tibble: 6 × 5
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## <dbl> <dbl> <dbl> <dbl> <fctr>
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
Conhecer o algoritmo de Árvore de Decisão e cirar um modelo que possa prever a espécie de Íris a flor a partir das medidas da pétala e sépala.
O primeiro passo é dividir nosso dataset em sets de treinamento e de teste. Para que as análises daqui pra frente sejam feitas somente no conjunto de treinamento e não influenciem as medidas feitas sobre o testset.
set.seed(654)
train_idx <- sample(nrow(iris), .75*nrow(iris))
train <- iris %>% slice(train_idx)
test <- iris %>% slice(-train_idx)
Ok, ao todo temos 4 preditores, vamos checar como os mesmos estão relacionados entre si.
train %>% ggpairs(aes(colour=Species), columns=1:4, lower="blank", diag="blank",
upper=list(continuous="points")) + theme_bw()
Dentro do espaço Preditor \(\Large{X}\) encotrar a melhor subdivisão que explique as classes, nesse caso as espécies de Íris.
\[ \Large{\mathcal{F}} : X \rightarrow Y \]
tree.fit <- rpart("Species ~ Petal.Length + Petal.Width", train, control=rpart.control(cp=0, minbucket=1))
fancyRpartPlot(tree.fit, sub="")
limiares <- rpart.subrules.table(tree.fit) %>%
as_tibble() %>%
select(Variable, Less) %>%
filter(!is.na(Less)) %>%
mutate(Less=as.numeric(as.character(Less)))
limiares_pl <- limiares %>% filter(Variable=="Petal.Length") %>% .$Less
limiares_pw <- limiares %>% filter(Variable=="Petal.Width") %>% .$Less
x_axis <- train$Petal.Length %>% range()
y_axis <- train$Petal.Width %>% range()
n_points <- train %>% nrow()
train %>% ggplot(aes(x=Petal.Length, y=Petal.Width, colour=Species)) +
geom_point() + theme_bw() + scale_x_continuous(breaks=seq(10,1)) +
geom_line(aes(x=rep(limiares_pl[1], n_points), y=seq(y_axis[1], y_axis[2], length.out=n_points)), colour="black") +
geom_line(aes(x=seq(limiares_pl[2], x_axis[2], length.out=n_points), y=rep(limiares_pw[1], n_points)), colour="black") +
geom_line(aes(x=rep(limiares_pl[2], n_points), y=seq(y_axis[1], y_axis[2], length.out=n_points)), colour="black") +
geom_line(aes(x=seq(limiares_pl[2], x_axis[2], length.out=n_points), y=rep(limiares_pw[1], n_points)), colour="black") +
geom_line(aes(x=rep(limiares_pl[3], n_points), y=seq(y_axis[1], limiares_pw[1], length.out=n_points)), colour="black") +
geom_line(aes(x=seq(limiares_pl[3], x_axis[2], length.out=n_points), y=rep(limiares_pw[2], n_points)), colour="black")
Cada região dessa representa uma folha. O espaço preditor foi segmentado aqui, porém se pusermos um exemplo em 3D como o de abaixo, vemos a estratificação.
Um exemplo aleatório:
Outro detalhe importante sobre a divisão do espaço de preditores é o que não acontece, como abaixo:
– Ideia Básica:
Como escolher o melhor preditor feature no nó da árvore?
Taxa do erro? -> Não converge bem, não é sensitiva ao crescimento da árvore.
Índice Gini -> Medida de pureza do nó.
\[G({s_i}) = \sum_{k=1}^{K}{\hat{p}_{{s_i}k}(1-\hat{p}_{{s_i}k})}\]
\[E({s_i}) = -\sum_{k=1}^K{\hat{p}_{{s_i}k}\log{{\hat{p}}_{{s_i}k}}}\]
Comportamento da função em relação a probabilidade das classes:
f_gini <- function(p){ p*(1-p) + (1-p)*(1-(1-p)) }
f_entr <- function(p){ ifelse(p%in%c(0,1), 0,
- (p*log(p, base=2) + (1-p)*log((1-p), base=2)))}
ps <- seq(0, 1, length.out=100)
y_gini <- sapply(ps, f_gini)
y_entr <- sapply(ps, f_entr)
ggplot(tibble(probs=ps)) +
geom_line(aes(x=probs, y=y_gini, colour="Gini")) +
geom_line(aes(x=probs, y=y_entr, colour="Entropia")) +
theme_bw()
1 - Medição de entropia no Nó Raiz, no caso usando somente os preditores: Petal Length e Petal Width.
O espaço Preditor é a região mostrada abaixo, no caso \(S_1\).
S1 <- train %>% select(Petal.Length, Petal.Width, Species)
ggplot(S1, aes(x=Petal.Length, Petal.Width)) +
geom_point(aes(colour=Species)) +
theme_bw()
Para calcular a entropia, precisamos somente da probabilidade de cada classe.
\[E({s_i}) = -\sum_{k=1}^K{\hat{p}_{{s_i}k}\log{{\hat{p}}_{{s_i}k}}}\]
S1 %>% group_by(Species) %>%
summarise(quantidade=n()) %>%
mutate(prob=quantidade/sum(quantidade)) %>%
mutate(prob=round(prob, 3))
## # A tibble: 3 × 3
## Species quantidade prob
## <fctr> <int> <dbl>
## 1 setosa 38 0.339
## 2 versicolor 38 0.339
## 3 virginica 36 0.321
Portanto a entropia daquele nó, pode ser calculado como:
E_S1 <- - (0.339*log(0.339, base=2) + 0.339*log(0.339, base=2) + 0.321*log(0.321, base=2))
print(E_S1)
## [1] 1.584349
2 - Critério para divisão:
Ganho de Informação: Soma ponderada da entropia dos nós filhos ser menor do que a do nó pai.
\[\Delta{E} = p(S_{1})E(S_{1}) - \Big[p(S_{11})E(S_{11}) + p(S_{12})E(S_{12})\Big]\]
No nosso caso temos duas variáveis contínuas. Para cada atributo, são testados todos os valores continuos do seu range. O limiar que obter o maior ganho de informação é o selecionado.
fs_entropy <- function(S){
if( nrow(S)==0){
return( 0)
}
list_p <- S %>%
group_by(Species) %>%
summarise(quantidade=n()) %>%
mutate(prob=quantidade/sum(quantidade)) %>%
.$prob
list_p <- list_p[list_p > 0 | list_p < 1]
entropia <- list_p %>%
sapply(X=., FUN=function(p){ -p*log(p, base=2) }) %>%
sum()
return(entropia)
}
delta_E <- tibble(attr=character(), thrs=numeric(), gain=numeric())
E_S1 <- fs_entropy(S1)
for( A in colnames(S1)[1:2]){
range <- S1 %>%
select(eval(parse(text=A))) %>%
range()
range_seq <- seq(range[1], range[2], 0.01)
for( i in range_seq){
S11 <- S1 %>% filter(get(A, pos=.) >= i)
S12 <- S1 %>% filter(!get(A, pos=.) >= i)
E_S11 <- fs_entropy(S11)
E_S12 <- fs_entropy(S12)
dE <- E_S1 - (1/nrow(S1))*(nrow(S11)*E_S11 + nrow(S12)*E_S12)
delta_E <- rbind.data.frame(delta_E, tibble(A, i, dE))
}
}
ggplot(delta_E, aes(x=i, y=dE)) +
geom_line(aes(colour=A)) +
facet_wrap(~A, nrow=1, scales="free_x") +
theme_bw()
Aqui vemos que existe um empate e os dois atributos atingem o valor máximo de ganho. Em petal.length esse valor de topo se repete dentro de uma faixa de valor.
delta_E %>% filter(dE==max(dE)) %>% filter(A=="Petal.Length") %>% select(i) %>% range()
## [1] 1.91 3.00
O valor médio dessa faixa, é exatamente o valor de corte escolhido na nossa primeira árvore. (O número aparente no diagrama árvore é arrendado por isso não bate examente.)
(1.91+3)/2
## [1] 2.455
3 - Repetir o passo 1 e 2.
Com os Nós filhos gerados da etapa anterior, repetir o processo de divisão.
\[\sum_{m=1}^{|T|}\sum_{x_i\in{S_m}}(y_i - \hat{y}_{R_m})^2 + \alpha|T|\]
Podando a árvore modelada anteriormente.
tree.fit.pr <- prune(tree.fit, cp=0.1)
fancyRpartPlot(tree.fit.pr, sub="")
Acertividade no trainset da árvore inteira e podada.
train.pred.1 <- predict(tree.fit, train, type="class")
## acertividade no trainset
mean(train.pred.1==train$Species) %>% round(3) %>% `*`(., 100)
## [1] 99.1
train.pred.2 <- predict(tree.fit.pr, train, type="class")
## acertividade no trainset
mean(train.pred.2==train$Species) %>% round(3) %>% `*`(., 100)
## [1] 95.5
Acertividade no testset da árvore na íntegra e podada.
test.pred.1 <- predict(tree.fit, test, type="class")
## acertividade no testset
mean(test.pred.1==test$Species) %>% round(3) %>% `*`(., 100)
## [1] 94.7
test.pred.2 <- predict(tree.fit.pr, test, type="class")
## acertividade no testset
mean(test.pred.2==test$Species) %>% round(3) %>% `*`(., 100)
## [1] 94.7
Construções de modelos bem mais poderosos, usando árvores.
Bagging
Boosting
Random Forest